rm(list=ls())

# Author: ScorpioGirl
# Date: 2019-12-07
# Aim: Download miRNA isoform expression from TCGA database
#install.packages("pacman")

##加载所需要的包以及环境配置
library(pacman)
p_load(TCGAbiolinks)
p_load(stringi)
packageVersion("TCGAbiolinks")

rm(list = ls())
options(stringsAsFactors = F)
setwd("C:\\Users\\lenovo\\Desktop\\data")

##数据查询与下载

# 可以通过TCGAbiolinks:::getGDCprojects()$project_id得到每个项目的ID
TCGAbiolinks:::getGDCprojects()$project_id
#project_id <- "TCGA-LUAD"
#project_id <- "TCGA-LUSC"
# project_id <- "TCGA-GBM"
#project_id <- "TCGA-LIHC"
#project_id <- "TCGA-KIRC"
project_id <- "TCGA-COAD"

# 可以通过TCGAbiolinks:::getProjectSummary(project_id)得到可以下载的数据类型
TCGAbiolinks:::getProjectSummary(project_id)
data_category <- "Transcriptome Profiling"

# data.type
data_type <- "Isoform Expression Quantification"

# 查询数据
query <- GDCquery(project = project_id,
data.category = data_category,
data.type = data_type)

# 下载
GDCdownload(query)


# 整合，这里得到的dataAssy是一个数据框，并且每一个样本的每一个miRNA一行，后面需要处理一下
dataAssy <- GDCprepare(query = query,
summarizedExperiment=F)

##得到有多少个样本

sample_id <- as.data.frame(table(dataAssy$barcode))
sample_id

##得到有多少个成熟体miRNA
miRNA_id <- as.data.frame(table(dataAssy$miRNA_region))
miRNA_id

##生成一个有xx个样本与xx个成熟体的miRNA矩阵
miRNA_matrue_RPM <- matrix(NA,ncol = nrow(sample_id),nrow = nrow(miRNA_id))
colnames(miRNA_matrue_RPM) <- sample_id$Var1
rownames(miRNA_matrue_RPM) <- as.character(miRNA_id$Var1)

##处理数据得到RPM值
#未表达的miRNA这里处理成了0，后续根据需要可以进行过滤

for(i in 1:nrow(sample_id)){
 temp1 <- dataAssy[which(dataAssy$barcode==as.character(sample_id[i,1])),]
 
 for(j in 1:nrow(miRNA_id)){
     loc <- which(temp1$miRNA_region==as.character(miRNA_id[j,1]))
     if(length(loc)>0){
       miRNA_matrue_RPM[j,i] <- sum(temp1[loc,4])
     }
     else{
       miRNA_matrue_RPM[j,i] <- 0
           }
     }
 print(i)
}

write.table(miRNA_matrue_RPM,"miRNA_matrue_RPM.txt",sep="\t",row.names=T,quote=F)

##需要进行修改
miRNA_matrue_RPM1 <- miRNA_matrue_RPM[1:2102,]


#对名字进行处理
#name <- unlist(strsplit(rownames(miRNA_matrue_RPM1),split="mature,",fixed=T))[seq(2,2166*2,by=2)]
write.table(miRNA_matrue_RPM1,"miRNA_matrue_RPM1.txt",sep="\t",row.names=T,quote=F)

name <- substr(rownames(miRNA_matrue_RPM1),8,nchar(rownames(miRNA_matrue_RPM1)))


# 使用R包miRBaseVersions.db进行ID转换
#miRNA_matrue_RPM1 <- read.table("miRNA_matrue_RPM1.txt")
library(miRBaseVersions.db)
items <- select(miRBaseVersions.db, 
                 keys = name, 
                 keytype = "MIMAT", 
                 columns = c("ACCESSION", "NAME", "VERSION"))

id_name <- items[items$VERSION == 21.0, c("ACCESSION","NAME")]

#未经尝试
#id_name <- read.table("id_name.txt")

miRNA_matrue_RPM2 <- cbind(id_name,miRNA_matrue_RPM1)

#可获得成熟体表达谱

write.table(miRNA_matrue_RPM2,"miRNA_matrue_RPM2.txt",sep="\t",row.names=F,quote=F)

